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Abstract. The line emission from a growing number of Herbig-Haro jets can be observed and resolved at angular 
distances smaller than a few arcseconds from the central source. The interpretation of this emission is problematic, 
since the simplest model of a cooling jet cannot sustain it. It has been suggested that what one actually observes 
are shocked regions with a filling factor of ~ 1%. In this framework, up to now, comparisons with observations 
have been based on stationary shock models. Here we introduce for the first time the self-consistent dynamics of 
such shocks and we show that considering their properties at different times, i.e. locations, we can reproduce ob- 
servational data of the DG Tau microjet. In particular, we can interpret the spatial behavior of the [SII]6716/6731 
and [Nil] /[OI] 6583/6300 line intensity ratios adopting a set of physical parameters that yield values of mass loss 
rates and magnetic fields consistent with previous estimates. We also obtain the values of the mean ionization 
fraction and electron density along the jet and compare these values with those derived from observations using 
the sulfur doublet to constrain the electron density. 
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1. Introduction 

Herbig-Haro jets can be observed at very high spatial res- 
olution by HST and with the next generation of ground 
based optical telescopes such as VLTI/ AMBER the angu- 
lar resolution capabilities will be boosted from the actual 
fraction of an arcsecond up to the milliarcsecond range 
(Bacciotti 2004). It will soon be possible to look into the 
very first part of the jet as it emerges from the accretion 
disk or from the reflection nebula and resolve many jets 
in their radial extent. It will thus be feasible to perform a 
comparison of direct observations with the predictions of 
the acceleration models. Some jets (sometimes called 'mi- 
cro jets'), e.g. HH 30 (Bacciotti et al. 1999) or DG Tau 
(Lavalley et al. 1997, Lavalley-Fouquet et al. 2000 (L- 
FCD2000), Bacciotti et al. 2000), are particularly good 
candidates for a high resolution study of the evolution of 
the physical parameters along the initial fraction of the 
jet, i.e. up to ~ 10 16 cm where the conditions to meet for 
the line emission mechanisms at work are the most severe. 

The jet of DG Tau has been observed close to its 
base at CFHT at high angular resolution (~ 0".5) by 
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L-FCD2000 who showed the behavior of [SII]6716/6731 
and [Nil] /[OI] 6583/6300 line ratios along the jet for the 
high (HV), intermediate (IV) and low (LV) velocity com- 
ponents. Our calculations will address these data in mod- 
eling the jet line emission. 

These observations typically show that the behavior of 
temperature, ionization and density along the jet is incom- 
patible with a freely cooling jet. Various heating processes 
have been proposed in the literature, such as ambipolar 
diffusion (Garcia et al. 2001a, b), photoionization by soft 
X-rays from the TTauri star (e.g. Shang et al. 2002) and 
mechanical heating (O'Brien et al. 2003, Shang et al. 2002, 
for X-wind jets). These estimates where carried out for 
steady-state jet models and pointed out that mechanical 
heating was the most effective in reproducing the obser- 
vations. The idea of tapping a small fraction of the jets 
kinetic energy to convert into heat is certainly appealing; 
however there is no physical explanation for if and how 
this process could work in YSO jets. Velocity fluctuations 
may possibly steepen into shocks and dissipate their en- 
ergy heating the gas, but that radiative losses come into 
play and act against the heating process. 

An alternative explanation by L-FCD2000 to inter- 
pret the line ratios of DG Tau (see also Hartigan (2004)), 
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is that one observes several post-shock regions of high 
excitation with a filling factor of <~1% ('shocking jet'). 
L-FCD2000 found that the line ratios of DG Tau and 
other HH objects (compilation by Raga et al. 1996) might 
be interpreted as series of shocks arrayed along the jet 
with varying shock velocities as one moves along the ra- 
dius away from the jet axis, typically: 70-100 km s _1 for 
the HV component, 50-60 km s _1 for IV and pre-shock 
densities that decrease from the star as oc r~ 2 starting 
from no = 10 5 cm~ 3 . In their analysis, they considered 
post-shock parameters consistent with planar, stationary 
shocks models (e.g. Hartigan et al. 1994). 

Raga et al. (2001) modeled the morphological and dy- 
namical properties of the DG Tau jet by means of 3D 
numerical simulations, assuming a precessing jet with a 
velocity that varied sinusoidally in time. In the present 
paper we restrict ourselves to the interpretation of the 
line intensity ratio behavior along the jet, abandoning the 
assumption of stationary shocks. We follow the dynamical 
evolution of an initial perturbation as it steepens into a 
(radiative) shock traveling along the jet, and derive the 
post-shock physical parameters consistently (Massaglia et 
al. 2005a). From these parameters we construct the emis- 
sion line ratios to be compared with observations. In this 
first analysis we consider the evolution of a single shock, 
neglecting the possible interaction with other shocks. 

The plan of the paper is the following: in Section 2 
we outline the initial conditions and the computational 
scheme adopted; in Section 3 we examine the shock evo- 
lution, while in Section 4 we discuss the results and make 
comparisons with observations. Our conclusions are drawn 
in Section 5. 

2. The model 

2.1. Basic equations 

We restrict our attention to one-dimensional MHD planar 
flow. The fluid is described in terms of its density p, veloc- 
ity u, thermal pressure p and (transverse) magnetic field 
By . Conservation of mass, momentum, magnetic field and 
total energy readily follows: 

dp d{pu) = Q 
dt dx 




+ P+ ■ 



= 



dB, 
dt 



v + d{B y u) = Q 



dx 



dE d 
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where E = p/{T-l)+pu 2 /2+B 2 /2 is the total energy den- 
sity (we use T = 5/3) and £(T, /„) represents the energy 
loss term (energy per unit volume per unit time) which 



depends on the temperature T and, as explained below, 
on the number fraction of neutral hydrogen atoms, /„. 
The loss term accounts for energy lost in lines and in the 
ionization and recombination processes. Line emissions in- 
clude contributions from nine elements whose abundances 
have been assumed to be solar (Ansers & Grevesse 1989): 
H and He resonance lines, the 13 strongest forbidden lines 
of C, N, O, S, Si, Fe and Mg. He is neutral whereas met- 
als are singly ionized and the abundance of C is 10% of 
the solar one. The ionization states of N and O are fixed 
to that of H by charge transfer. 

To compare between the observed and com- 
puted line ratios, we have computed the popula- 
tions of the atomic levels for the forbidden transi- 
tions [SII]AA6716, 6731, [NII]A6583 and [OI]A6300, solv- 
ing, according to Osterbrock (1974), the excitation - de- 
excitation equilibrium equations for five energy levels. 

An additional evolutionary equation is solved for the 
neutral fraction / n , 



^+^=n c [- Ci / n + Cr (l-/ n )] 



(5) 



where n c is the electron density, whereas c; and c r are, 
respectively, the ionization and recombination rate coeffi- 
cients (Rossi et al. 1997). In this framework we have 



"h(1 - /n) + Zrin 



and 



p = n H (2 - / n + Y + 2Z)k B T 



(6) 



(7) 



where tih is the total hydrogen density, Y and Z (= 0.001) 
are the Helium and metal abundances by number respec- 
tively. We will also define the ionization fraction as 



f _ n c 

Ji — 

n H 



(8) 



2.2. The initial perturbation 



A nonuniform pre-shock density that decreases away from 
the star is a reasonable ingredient when dealing with an 
expanding jet, as suggested by observations. Moreover, 
shocks that propagate into a stratified medium tend to 
increase their amplitude when they find a decreasing 
pre-shock density. Therefore, we consider the following 
nonuniform pre-shock density: 



p {x) = po 



xi + x l 



(9) 



where x is the spatial coordinate along the jet axis. Thus, 
for a; > iq, we have a conical decrease of density, i.e. a 
conical expansion of the jet, while for x <x n the density 
decreases parabolically, po(x) ~ po(l — x 2 /xq). x sets the 
initial steepness of the density function and this affects 
the shock evolution even at larger distances. 

A generic initial perturbation imposed above the mean 
flow typically evolves forming two shocks: the forward and 
reverse shocks (Hartigan & Raymond 1993, Massaglia et 
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al. 2005a). The energy content of the perturbation there- 
fore splits into two radiating shocks that propagate, de- 
creasing their strengths along the way. Thus one cannot 
reach and maintain the compression factors to explain ob- 
servations, as happened in Massaglia et al. (2005a), un- 
less assuming an initial perturbation velocity amplitude 
of the order of the jet mean flow speed. Thus it would 
be more desirable to work with a single forward shock. 
To this purpose, we consider a disturbance that maintains 
the Riemann invariant J_ constant (Zeldovich & Raizer 
1966): 



J- 



dp 
pc 



dp du 

— = constant 

du pc 



After differentiation and some algebra, one finds 
7-1 



= (u- U ) +p : 



(10) 



(11) 
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We set Uq(x) = 0, implying that we carry out our calcula- 
tions in the reference frame of the mean flow. We will have 
to transform our data to the laboratory frame when com- 
paring results with observations. We prescribe the velocity 
perturbation as 



u(x) 







-{x — xi) 2 + 2a(x - 



x\)] if 2a + x\ > x > X\ 
otherwise 




where x\ is the initial coordinate of the perturbation and 
a is its half-width (see Fig.^ solid line). Different choices 
of the initial perturbation shape are not crucial for the 
shock formation and evolution. 

Notice that, strictly speaking, J_ in Eq. is no 
longer an hydrodynamic invariant when magnetic fields 
are present. In our case, however, the initial magnetic 
pressure is typically small and the resulting perturbation 
still produces a strong forward shock, as shown in Fig. ^ 
In our calculations we have assumed x\ = 10 14 cm and 
cr = 2 x 10 13 cm. 

The boundary conditions assume free outflow at x = 
and x = L. The extent of the computational domain, 
L, has been chosen sufficiently large to follow the shock 
evolution for t ~ 15 ys and avoid spurious interactions 
with the boundaries. For this reason we adopt L = 4.5 x 
10 15 cm. 



2.3. Method of solution 

To solve the MHD equations we have employed the hy- 
drodynamical code PLUTO (Mignone, Massaglia & Bodo, 
2004). 

Equations Q— Q, together with Eq. JSJ are solved us- 
ing a conservative, second-order accurate total variation 
diminishing (TVD) scheme. Piecewise linear interpolation 
with limited slopes is used to ensure monotonicity inside 
each computational zone. Second order accuracy in time 
is achieved via a MUSCL-Hancock predictor step (van 
Leer, 1985). A linear, Roe- type, Riemann solver is used 
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Fig. 1. Density (upper panel) and velocity (lower panel) 
vs distance in the reference frame at rest with the mean 
flow (abscissa x m f); solid lines show the initial perturba- 
tion, the dashed lines the evolution at t — 5 ys and dot- 
dashed lines at t = 10 ys 



to compute the inter-cell fluxes needed in the conserva- 
tive update. The conservative formulation is essential for 
a correct description of the shock dynamics (LeVeque et 
al., 1998). Source terms describing cooling, ionization and 
recombination processes are treated using operator split- 
ting. 

In the simulations presented below, the region of in- 
terest is confined mainly to the post-shock flow, where 
most of the emission takes place. Since this region is much 
smaller in size (< 7 x 10 13 cm) than the domain and 
the shock front is not stationary, a static uniform grid 
would demand increasingly high resolution in order to ad- 
equately resolve the post-shock structure. For this rea- 
son, we use an Adaptive Mesh Refinement (AMR) tech- 
nique (described in Massaglia et al. 2005b) when solving 
Eqs. 0-®, thus reaching an effective resolution of about 
360, 000 grid points. 

This allows to considerably speed up the computation 
with sufficient resolution to accurately describe the dy- 
namics and emission processes. 
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3. Results and discussion 

We have followed the temporal evolution of the perturba- 
tion as it steepened into a shock traveling along the com- 
putational domain (see Fig. We have been able, with 
the AMR technique employed, to resolve the post-shock 
region at all times and in Fig. [21 we show the post-shock 
variables versus the spatial coordinate. Our model cal- 
culations depend on the following parameters: the initial 
perturbation velocity amplitude Uo, the pre-shock mag- 
netic field transverse to the flow Bq, the pre-shock den- 
sity parameters Xq and no and the mean flow speed Uq, 
to transform the results back to the observer's reference 
frame. Test calculations have shown that, among these 
parameters, uq and the pre-shock density parameters xq 
and no are the most critical for the final results. We have 
fixed the (isothermal) pre-shock temperature T — 1,000 
K (Raga et al. 2001) and ionization fraction f l = Z due to 
metals only. However results are weakly sensitive to these 
values, for the range of parameters considered. In Fig.[2]we 
plot the post-shock temperature T, electron density (n e ), 
ionization fraction (/i), magnetic field B y in units of Bq, 
[SII] emissivity (in units of the maximum value, rightmost 
vertical scale) and velocity u, in the reference frame of the 
mean flow. The horizontal scale has an arbitrary origin 
and tells us the size of the post-shock region. The values 
of the parameters are: u = 70 km s~ 4 , B = 100 /zG, 
x = 0".l and ti = 5x 10 4 cm~ 3 . 

The post-shock quantities are shown after 10 ys of 
evolution, when the shock front has traveled for about 
6.5 x 10 15 cm from the source, having assumed a mean 
flow speed U = 150 km s -1 (see below). We note that 
the post-shock temperature attains a value > 10 5 K right 
behind the shock and decreases to below 10 4 K about 
7 x 10 13 cm from the shock front, where the [SII] emission 
takes place and the magnetic field reaches its maximum 
and subsequently drops to low values at ~ 10 14 cm. We see 
that both the (collisional) ionization fraction and the elec- 
tron density reach maximum values at about 7 x 10 13 cm 
behind the shock front. In the post-shock region the par- 
ticle velocity remains nearly constant to the original value 
of the perturbation amplitude and decreases starting at 
~2x 10 15 cm from the shock front. We are not plotting 
the post-shock speed in the frame where the shock is at 
rest. At earlier times, the general behavior would be al- 
most unchanged for each quantity exception made for the 
electron density which would present a maximum that is 
higher by about a factor of three due to the higher (im- 
posed) pre-shock density. Clearly, the opposite happens at 
later times in the evolution. 

In order to carry out a comparison with observa- 
tions we need to average the post-shock quantities at ev- 
ery evolutionary time point. We have done this following 
Hartigan et al. (1994), i.e. defining the [S'//]-weighted av- 
erage as follows 

jQ(x)e{[SII](x)}dx 
W> fc{[SII](x)}dx 
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Fig. 2. Spatial behavior, in the mean flow frame, of the 
physical quantities as they undergo the shock transition, 
from left to right, after lOys of evolution, corresponding to 
a traveled distance of the shock front of about 6.5x 10 15 cm 
along the jet axis. Labels are: f\ ionization fraction (solid 
line), n c electron density (dashed line), T temperature 
(dotted line), u velocity (dot-dashed line), B y magnetic 
field (in units of the pre-shock value B , dot-dot-dashed 
line) and [SII] emissivity (solid line) 



where Q is a physical quantity such as either electron 
density or ionization fraction or the line emissivity ratio 
of the sulfur doublet [57/]A6716/[5//]A6731. Concerning 
the line emissivity ratio [iV//]A6583/[O/]A6300, we have 
carried out the averaging procedure in the same way, but 
adopting the total emissivity e{[NII}(x)} + e{[OI](x)} as 
the weighting function. However, we noticed that the par- 
ticular choice of the weighting function has very little ef- 
fect on the results. Note also from Fig. that for /j and 
By the averaging procedure will yield, with a good ap- 
proximation, the maximum value of these quantities. 

In Fig. El we compare of the averaged post-shock quan- 
tities, obtained for the above set of parameters and assum- 
ing a mean outflow velocity Uq = 150 km s _1 , with the 
observations by L-FCD2000. The abscissa indicates the 
position of the post-shock region in the observer's frame, 
x = Uot + £(t), where t is the elapsed time and £(t) is the 
location of the emission region. We have chosen the high 
velocity (HV) jet data of the line ratios, ionization fraction 
and electron density from the paper of L-FCD2000, likely 
corresponding to the inner part of the jet about its longitu- 
dinal axis. After examining Fig. [21 we have estimated 
as the shock front position minus an interval d (= 2 x 10 14 
cm) that accounts for the shift of the emission region be- 
hind the shock front. The results are nearly insensitive to 
the choice of d. We recall that the directly observable data 
are the line intensity ratios (in Fig. [31 crosses are for the 
[SII] doublet and diamonds for the [Nil]/ [OI] ratio), while 
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Fig. 3. Averaged values of the post-shock electron density 
((n c ), dashed line), ionization fraction ((/;) solid line), 
ratios of [SII] (dot-dashed line) and [Nil]/ [01] (dot-dot- 
dashed line) along the jet for the reference parameters 
vs distance (in the laboratory frame): uq = 70 km s _1 , 
Bq = 100 //G, x = 0".l and n = 5 x 10 4 cm -3 (see also 
Fig.|2J). Here we have set Uq = 150 km s _1 



the ionization fraction and electron density have been ob- 
tained by L-FCD2000 adopting the diagnostic technique 
described in Bacciotti et al (1995). From Fig. [3|we note 
that the [Nil]/ [OI] and [SII] ratios fit the data remark- 
ably well, while the electron density and the ionization 
fraction exceed the ones derived from observations at low 
values of the spatial coordinate x. In fact, Bacciotti et 
al. (1995) obtain n e from the sulfur doublet line intensity 
ratio, which saturates for electron densities higher than 
the critical density n c < 10 4 cm~ 3 . This leads to an un- 
derestimation of the electron density and, consequently, 
the ionization fraction. Note that in L-FCD2000 the first 
couple of data points of the electron density and the first 
one of the ionization fraction are lower limits. Instead, 
this method is appropriate for lower values of the electron 
density, thus the agreement in Fig. [3] at larger distances. 

Several questions arise: which sets of parameters better 
represent the observations? What are the crucial parame- 
ters? 

Variations of the advection velocity Uq have no phys- 
ical impact on the shock evolution but concern only a 
dilation of x axis. We varied Uq by ± 40% and noted that 
the agreement remained reasonably good. Thus values of 
Uq within these two limits still interpret observations rea- 
sonably well. 



Table 1. Behavior of the line ratio curves with the in- 
crease of the parameter in the leftmost column 



Parameter 


[Nil] 6583/ [OI] 6300 


[SII]6716/[SII]6731 


no 


I 


I 


Uo 


T 


I 


Xo 


I 


I 


B Q 


T 


T 



We are left with four main parameters: no, xq, Bq 
and uo- Due to the difficulty in carrying out a full ex- 
ploration of the parameter space, we have synthesized the 
'goodness' of the agreement between model curve and ob- 
servational points of the line ratios defining a variance 
a 2 = J2i (Oi — Ei) 2 , where Oi are the data and Ei the 
model points. In Table 1 we show the qualitative behav- 
ior of the [SII] and [NII]/[OI] ratio model curves when 
a parameter on the left column is increased, leaving the 
other ones unchanged with their reference values of Fig.|3J 
these curves are displaced downward when no and Xq are 
raised, changing in shape as well, while the curves shift 
upward when raising Bq. We have run several models (ac- 
cording to Table 1) varying the parameters no, -Bo, uq and 
Xq around the reference set and calculated the variance a 2 
of the corresponding fits. Since for most of the spatial co- 
ordinate x we are above the critical density, the variance 
of the fits of the [SII] doublet line ratio is not very sen- 
sitive to the changes in these parameters. However, this 
is not the case for the [NII]/[OI] curves: we could obtain 
fits for [NII]/[OI] with a variance close to the minimum 
by varying no, xq and Bq by about 30% and uq by about 
10% plus or minus with respect to the reference values. 

As far as the magnetic field evolution is concerned, di- 
rect comparisons with previous 2D models (e.g., Garcia 
et al. 2001a) and 3D simulations (e.g., Cerqueira & de 
Gouveia dal Pino 2001) are difficult because we only in- 
clude the transverse component of the magnetic field. The 
spatial evolution of the mean field in our models typi- 
cally shows a monotonic decrease from about 30£?o down 
to about 5£>o- If we consider the behavior of B^ from 
Garcia et al. (2001a), our values of the magnetic field, 
with Bq = 100/iG, remain a factor of ~ 2 — 3 below the 
curve of their model C. 

4. Conclusions 

With the goal to interpret the line-emission within the 
first 5" observed in some Herbig-Haro jets, we have con- 
sidered the evolution of a strong planar perturbation in 
a stratified, magnetized medium and in the presence of 
radiative losses. We have examined the temporal evolu- 
tion of this perturbation as it steepened into a shock and 
focused our attention to the structure of the post-shock re- 
gion. Having set physical parameters consistent with the 
environment of stellar jets in their inner portions, closer 
to the young star, we have derived averaged line intensity 
ratios that could be compared with observations. We have 
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adopted observations of line ratios of the DG Tau jet by 
L-FCD2000, as an example, and found that even in this 
extremely simple model, one can interpret observations 
with a constrained set of parameters reasonably well. 

In this model we have a priori prescribed the longitu- 
dinal density profile upstream of the perturbation, assum- 
ing that a real jet, in the same physical conditions, would 
expand decreasing its density along the way. Whether an 
actual jet behaves in a similar fashion remains to be seen 
and, to this purpose, 2D MHD simulations of radiative jets 
are needed. The reference set of parameters that yield a 
'good' agreement with observations are: uq = 70 km s _1 , 
B a = 100 /uG, x = 0".l, n = 5 x 10 4 cm~ 3 with a mean 
flow speed Uq = 150 km s _1 . Assuming an initial radius of 
3 x 10 14 cm (Raga et al. 2001), these values are consistent 
with a mass loss rate M w 5 x 10~ 8 M Q ys _1 , in agree- 
ment with the estimates of L-FCD2000. The presence of a 
perturbation raises this value to M « 10~ 6 M Q ys _1 for 
short periods. 

One might ask whether these shocks may also be re- 
sponsible for the emission knots observed in some HH jets 
(e.g. HH34, HH111) at larger distances from the source 
(< 45") (c.f. a review by Raga, Beck & Riera 2004, Micono 
et al. 1998a, b). Due to the strong emissivity of the shock 
compressed medium, we believe that these shocks would 
hardly be visible at very large distances, as observed. 

Thus, these calculations successfully reproduce line- 
ratio observations and thus strongly support the hypothe- 
sis of L-FCD2000 of the DG Tau jet as a 'shocking jet', i.e. 
that one actually observes not a continuous emitting jet 
but just the gas parcels that have undergone compression, 
heating and ionization in shocks. 
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